#ifndef  __CPARTSYS_H__     // Begin 'include-guard' block.
#define  __CPARTSYS_H__     // (See #endif at the very end of this file).

//==============================================================================
//3456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_
//                               (why all the numbers? I prefer 80-column text)
//
//  CPartSys.h  //  Master Particle System class.
//              One variable of this type holds an entire particle system.
//              SEE bottom of file for detailed description..
//==============================================================================

// multi-platform includes for the GLUT/OpenGL libraries:-----------------------
#if defined(_WIN32)
#include <GL\glew.h>        // OpenGL 'extensions wrangler' for easy access.
#include <GL\freeglut.h>// open-source GLUT functions for windows, keybd, mouse
/*---*OR*--- use:
// #include <GL/glut.h> // for 100% backwards-compatibility with Mark Kilgard's
                        // SGI-copyrighted GLUT library. Careful! this excludes
                        // freeglut's callbacks for mouse wheel, etc.
*/
/*
// freeglut.h already contains these other includes that were necessary for
// using Kilgard's original GLUT library:
#include <GL/GL.h>          // OpenGL 3D graphics rendering functions
#include <GL/Glu.h>         // OpenGL utilties functions
#include <GL/glext.h>       // #defines for all OpenGL features on your machine.
*/
#elif defined(__APPLE__)
#include <GLUT/glut.h>      // GLUT functions for windows, keybd, mouse
#include <OpenGL/GL.h>      // OpenGL 3D graphics rendering functions
#include <OpenGL/Glu.h>     // OpenGL utilties functions
#include <OpenGL/GLext.h>   // #defines for all OpenGL features on your machine.
// UNIX, Linux of all flavors;
#else
#include <unistd.h>         // defines symbolic constants, POSIX versions, etc.
#include <GL/glew.h>        // OpenGL 'extensions wrangler' for 'mesa' lib usage
#include <GL/freeglut.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glext.h>
#endif
//-----------------------------end multi-platform includes for GLUT/OpenGL-----

// C++ math, I/0 libraries (defined in the 'std' namespace)====================
#include <cmath>            // For sqrt, pow, etc.
#include <iostream>         // for cin, cout...
using namespace std;        // to access iostream's namespace

#include "CPart.h"              // CPart class; one individual particle.
#include "CForcer.h"            // CForcer class: one force-applying object
#include "CWall.h"              // CWall class: one constraint-applying object.

// Literals to select what kind of solver to use in CPartSys:
//------------------------------------------------------------------------------
// EXPLICIT methods: GOOD!
//          ++ simple, easy to understand, fast, but
//          -- Requires tiny time-steps for stable stiff systems, because
//          -- Errors tend to 'add energy' to any dynamical system, driving
//              many systems to instability even with small time-steps.
#define NU_SOLV_EULER       0       // Euler integration: forward,explicit,...
#define NU_SOLV_MIDPOINT    1       // Midpoint Method (see Pixar Tutorial)
#define NU_SOLV_ADAMS_BASH  2       // Adams-Bashforth Explicit Integrator
#define NU_SOLV_RUNGEKUTTA  3       // Arbitrary degree, set by 'solvDegree'

// IMPLICIT methods:  BETTER!
//          ++Permits larger time-steps for stiff systems, but
//          --More complicated, slower, less intuitively obvious,
//          ++Errors tend to 'remove energy' (ghost friction; 'damping') that
//              aids stability even for large time-steps.
//          --requires root-finding (iterative: often no analytical soln exists)
#define NU_SOLV_IMPLICIT    4       // Iterative implicit solver;'back-winding'

// OR SEMI-IMPLICIT METHODS: BEST?
//          --Permits larger time-steps for stiff systems,
//          ++Simpler, easier-to-understand than Implicit methods
//          ++Errors tend to 'remove energy) (ghost friction; 'damping') that
//              aids stability even for large time-steps.
//          ++ DOES NOT require the root-finding of implicit methods,
#define NU_SOLV_VERLET      5       // Verlet semi-implicit integrator;
#define NU_SOLV_VEL_VERLET  6       // 'Velocity-Verlet'semi-implicit integrator
#define NU_SOLV_MAX         7       // number of solver types available.

#define NU_EPSILON 10E-15           // tiny amount; a minimum vector length
                                    // to use to avoid 'divide-by-zero'

class CPartSys
//==============================================================================
// A complete particle system, consisting of a group of particles, a group of
// forces, a group of numerical integrators to find how they change with time.
{
public:
// Operate this Particle-System object:
//---------------------------------------
    CPartSys(void);                 // Default constructor/destructor
    ~CPartSys(void);

    // Initialize: call one of these fcns in main() (or in keyboard callbacks).
    void initExample1(void);        // Cube of particles outside a bubble.
    void initExample2(void);        // Cube, bubble, gravity and drag.
    void initExample3(void);        // Springy objects, bubble, gravity, drag

    // Run: call this in GLUT_registered myDisplay() callback,
    // or use myIdle() callback to control animation as I did here.
    void solver(void);              // Use the selected solver type to find the
                                    // next state for the particle system in s1;
                                    // then swap s1 <-> s0 and s1dot <->s0dot.
    void solverOLD(void);   // temporary--for removal of CMotives

    // Draw: call this in GLUT-registered myDisplay() callback,
    void drawMe(void);              // Render the current particle system state.


    //----------solver() Support-------------
    GLdouble timeStep;              // Timestep used by the solver(); init'd
                                    // in the CPartSys constructor.
    int solverType;                 // Type of numerical integration method used
                                    // to compute s1 from s0,s0dot, and movtive.
                                    //  see #define NU_SOLV_... statements above
    int solvDegree;                 // Used only for Runga-Kutta solver: sets
                                    // the number of Taylor Series terms for the
                                    // polynomial fit: 2==quadratic, 3==cubic,...
    void dotMaker(CPart *pDotDest, CPart *pSrc, CForcer *pF);
                                    // Find pDotDest, time-derivative of the
                                    // state vector pSrc.  The derivative is the
                                    // one thing we can find exactly; it is the
                                    // result of Newton's Law (F=ma) and all the
                                    // force-making objects in the list at pF.
                                    // (how many? CPartSys::forceCount )
    void applyAllForces(CPart *pS, CForcer *pF);
                                    // Complete the state pS; accum. the total
                                    // force applied to all the particles by
                                    // all the forces in the list at pF.
    void stateVecSwap(CPart **ppA, CPart **ppB);// swap state-vector contents by
                                                // swapping pointer values.
    //--------Animation Run/Stop------------
    bool runNotStep;                // Animation Control:single-step when false;
                                    //                  run animation when true;
                                    // used by GLUT callback functions myIdle(),
                                    // myDisplay() and myKeyboard() function to
                                    // respond to 'p' key and space-bar key.

    //-------State Vectors-----------------------!
    int partCount;                  // # of particles held in each state vector.
                                    // (if <=0, state vectors set to NULL)
    CPart *pS0, *pS0dot;            // s0: state vector holding CURRENT state of
                                    // particles, and its time derivative s0dot.
    CPart *pS1, *pS1dot;            // s1 state vector holding NEXT state, and
                                    // its time-derivative.
    CPart *pSM, *pSMdot;            // sM, sMdot: midpoint state vector and its
                                    // time-deriv.(see solver, midpoint method)
                                    // (NOTE: same size for all state-vectors).
    // OPTIONAL: as needed,
    //      create additional state-vectors to implement other solvers such as
    //      Verlet, Adams-Bashforth, etc.
    // !!CAREFUL!! Be sure to update stateVecSizer(), stateVecWipe(),
    //             CPartSys(), and ~CPartSys() for the new state-vectors;
    //  NOTE: all state-vectors are the same size (CPartSys::partCount).
    void stateVecSizer(int nuSize);// For each state vector, create a new array
                                    // (dynamic allocation) of particles. Discard
                                    // any previous particles (if any) and if
                                    // successful, set partCount to 'nuSize'.
                                    // (NOTE: same size for all state-vectors).
    void stateVecWipe(void);        // Discard all existing memory allocated for
                                    // state vectors in this CPartSys.

    void stateVecDraw(CPart *pS);   // Draw state-vector pS usingOpenGL calls.

    // Math operators (BETTER: use operator overloading!)
    void stateVecZero(CPart *pS);                   // set state to zero
    void stateVecAssign(CPart *pDest, CPart *pSrc); // copy src to dest
    void stateVecAplusBtimesC(CPart *pDest, CPart *pA, CPart *pB, double sc);
                                                    // dest = A + sc*B


    //-------Force-List Vector(s)----------------
    int forcerCount;                // # of forcer-making objects held in the
                                    // dyn. alloc'd list at pF0 (if needed, in
                                    // pF1, pFM, pF0dot,pF1dot,pFMdot as well).
    CForcer *pF0;                   // f0; forcer-vector-- dyn. alloc'd list of
                                    // all CURRENT force-making objects,
    // *pF0dot;                     //  and its time-derivative f0dot
                                    // (==NULL if none).
/*    // OPTIONAL: as needed,
    //      create these additional forcer-vectors to implement time-dependent-
    //      and particle-dependent (?and wall-dependent?) force-making objects.
    //  !!CAREFUL!!! Be sure to update forcerVecSizer(), forcerVecWipe(),
    //               CPartSys() and ~CPartSys() for the new forcer-vectors;
    //  NOTE: all forcer-vectors are the same size (CPartSys:: forcerCount).
    CForcer *pF1, pF1dot;           // f1; forcer-vector--dyn. alloc'd list of
                                    // all NEXT force-making objects, and its
                                    // time-derivative f1dot (==NULL if none)
    CForcer *pFM, *pFMdot;          // fM, fMdot; midpoint forcer-vector and its
                                    // time deriv. (see solver, midpoint method)
*/
    void forcerVecSizer(int nuSize);// Change the size of each dyn.-alloc'd list
                                    // of force-applying objects (CForcer obj.).
                                    // Discard any previous forcers (if any),
                                    // and if successful, set forcerCount to
                                    // nuSize.(same size for all forcer-vectors)
    void forcerVecWipe(void);       // Discard all existing memory allocated
                                    // for forcer-vectors in this CPartSys.
    void forcerVecDraw(CForcer *pF, CPart *pS);
                                    // Draw forcer-vector pF using OpenGL calls;
                                    // (needs corresponding state-vector at pS
                                    // because force-applying springs get their
                                    // positions in space from the particles
                                    // they connect together).
    //-------Wall-Vectors(s)---------------
    int wallCount;                  //# of constraint-making objects (CWall obj)
                                    // held in the dyn. alloc'd list at pC0.
    CWall *pC0;                     // c0; constraint-vector--dyn. alloc'd list
                                    // of all CURRENT constraint-making objects,
    //CWall *pC0dot;                  // and its time-derivative c0dot
                                    // (==NULL if none).
    // OPTIONAL: as needed,
    //      create these additional wall-vectors to implement time-dependent-0
    //      and particle-dependent (?and forcer-dependent?) wall-making objects.
    // !!CAREFUL!! Be sure to update wallVecSizer(), wallVecWipe(), CPartSys()
    //      and ~CPartSys() for the new forcer-vectors.
    // NOTE: same size for all wall-vectors)
    void wallVecSizer(int nuSiz);   // For each limit-vector, create (by dyn.
                                    // allocation) an new array of constraints
                                    // (CLimit objects). Discard any prev.
                                    // constraints (if any), and if successful,
                                    // set limitCount to nuSiz.
                                    // (NOTE: same size for all limit-vectors).
    void wallVecWipe(void);         // Discard all existing memory allocated
                                    // for wall-vectors in this CPartSys.
    void wallVecDraw(CWall *pC);    // Draw the wall-vector pC with OpenGL calls

};


//==============================================================================
//
//	OpenGL/ GLUT 'starter' code meant as a flexible testbed for any and all
//      kinds of 3D particle systems, organized as taught in Jack Tumblin's
//      EECS 351-2 Intermediate Graphics class at Northwestern University.
//
//  C++ classes are organized for:
//  --Multiple independent particle systems (CPartSys objects), and
//  --Selectable, state-space solvers: these perform integration & timesteps
//        in different ways, such as explicit (Euler or 'forward') methods,
//        midpoint method, Runga-Kutta, implicit, Verlet, Adams-Bashforth, etc
//
//  --Each particle system object holds:
//      CPart objects:-----------------
//      --Individual descriptions of each particle (CPart objects) that each
//        include the same set of parameters, such as: position, velocity, mass,
//        charge, mass, age, fatigue, fear, and any other behavior-affecting
//        attributes of the particle itself that you may wish to add.
//      --state-vector descriptions of groups of particles; each state vector
//        is a pointer-to-CPart variable in CPartSys that holds the address of
//        a dynamically allocated array of CPart objects). All state-vectors
//        have the same size (the same number of particles==CPartSys::partCount)
//        and each CPart object has the same internal organization, and thus the
//        same stride between CPart members in each state vector.
//      CForcer objects:----------------
//      --Individual descriptions of each force-applying object(CForcer objects)
//        to impose any desired external force-like efx on any particle, such as
//          --Earth's gravity (makes things fall 'down'),
//          --Planetary gravity (massive particles attract others),
//          --Viscous drag (force = -velocity*K_drag)
//          --Electrostatics (charged particle attraction/repulsion),
//          --Springs, each one connects just 2 particles (a network of springs
//          can make a shape; rope, cloth, jello), etc.
//      --'forcer-vector' descriptions of groups of force-applying objects; each
//        force-vector is a pointer-to-CForcer variable in CPartSys that holds
//        the address of a dynamically allocated array of CForcer objects (e.g.
//        pF0, etc).  All forcer-vectors have the same size (the same number of
//        CForcer objects == CPartSys::forcerCount).  Each CForcer object in a
//        forcer-vector has the same internal organization, and thus the same
//        stride between CForcer members in each forcer-vector.
//      CWall objects:------------------
//      --Individual descriptions of each constraint-applying object(CWall obj.)
//        to impose any desired absolute constraint on the position or any other
//        state parameter on any particle, such as:
//          --a 'wall' or hard surface that no particle can pass through,
//          --a 'stick' or 'link' that keeps a specified pair of particles at a
//              fixed separation distance (e.g. www.SodaPlay.com links)
//          --a hinge, axle, piston, sliding joint, or other limitation on the
//              movement of one or more particles,
//          --a 'motor' that sets the position, velocity, or some/any other
//              parameter of a particle or set of particles, etc.
//      --'wall-vector' descriptions of groups of constraint-applying objects;
//        each wall-vector is a pointer-to-CWall variable in CPartSys that holds
//        the address of a dynamically allocated array of CWall objects
//        (e.g. pC0). All wall-vectors have the same size (the same number of
//        CWall objects==CPartSys::wallCount).  Each CWall object in a
//        wall-vector has the same internal organization, and thus the same
//        stride between CWall members in each wall-vector.
//
//  USEFUL TRICKS:
//    1) Particles hold parameter arrays and #define statements.  Instead of
//      uniquely named individual variables (e.g. mass, xpos,ypos,zpos, etc)
//      I use a fixed-size array of doubles, and use #define statements to
//      choose index #s for each unique name.  (SEE CPart.h)
//    2) 'partType': switch sign to enable/disable
//    Dynamic memory allocation is slow, and template 'collections' of objects
//      can be complex and error prone. I don't actually create/destroy CPart
//      CForcer, and CWall objects while the particle system runs: instead,
//          --Use dyn. alloc'd arrays of CPart, CForcer, CWall objects;
//          --Use an integer 'type' member in CPart,CForcer, CWall class that
//              ==0 for unused, 'dead' objects, and use other #define'd values
//              that specify what kind of object (e.g. CPart::partType tells us
//              how to render the particle; CForcer::forceType indicates it is
//              gravity, viscosity, spring, etc., CWall::wallType indicates it
//              is a wall, a stick, a hinge, a pulley, etc.).
//          --To TEMPORARILY disable/enable a CPart, CForcer, CWall object, just
//              change the sign of the 'type' member, not the value.
//          --If you need to 'create' a new particle or force, just search
//            the pS0, pF0, pC0 vectors to find a 'dead' one to use.
//  3) 'Hysteresis' in creating new particles greatly reduces time spent on
//      malloc(), free() (new, delete) and 'bookeeping' for particles:
//   --If you need more particle, forcer, or wall objects, don't create new ones
//      one-at-a-time; instead DOUBLE the size of the vectors, copy the old
//      objects to the new vector, and set all the unused objects' 'type' member
//      to 0, or 'dead'/unused.
//   --If you think you have TOO MANY particle, forcer, or wall objects, so
//      many that it seems you should delete them,
//          a) resist the urge--memory is cheap and plentiful! or,
//          b) wait--run a counter and look for, say, 128 successive requests to
//              reduce unnecessary particles before you act, or
//          c) wait: don't make a particle-reduction request until the number of
//              unused particles > 4X the number of particles in use, and then
//              cut # of particles in half.  This policy of grow / shrink by
//              2X / 0.5X; don't shrink until count persists below 0.25Max)
//              will usually keep resizing events rare, even for very dynamic
//              particle systems.
//==============================================================================
/*
CRUDE TUTORIAL------------------------
Particle System works in a few huge steps repeated over and over:
    Initialize: define the contents of a state vector s0, which describes the
    'situation' or 'state' of ALL particles in a single vector. Define the
    contents of f0, a collection of force-making objects (CForcer) that
    push and pull on particles based on their current state s0. Define the
    contents of w0, a collection of constraint-making objects (CWall) that
    prevent particles from violating rules such as "Stay at or above the floor,
    and do not pass through solid shapes", or "Keep fixed separation distance
    between these two particles connected by a 'stick'".

    Repeat these steps
    1) DOTFINDER: Find current sDot from current s and current forces
     --Given a state vector s0, find its derivative sDot by:
        --finding the net forces on each particle by summing up all the
            forces that apply to it in its current state.
        --using Newton's law: force/mass = acceleration = sDot for velocity
        --using other behavior laws: current rate-of-change for mass, etc.
        stateDot = DOTFINDER(state, forceLaws)
    2) SOLVER: find NEW s from current s and current sDot and any other useful
        state variables, forces, etc.
    3) SWAP: replace current s with new s (makes sDot, net forces invalid)
        draw particles in this newly-updated state, go back to 1).

***Why use a FORCE ACCUMULATOR in CPart?***
The 'apply forces' function (within the 'dotMaker' fcn) can be written several
different ways, but one has particularly strong advantages.  For every particle,
we need to find the vector total of the forces from every force-making object
that affects that particle, suggesting two nested loops; we could either
--step through each particle, and find the total of all the forces on that
particle, then use Newton's law to find the missing s0dot values (e.g. the
time-rate-of-change of velocity).
--step through each FORCE, and find the force it applies to any/all particles.
This approach ('loop-thru-forces' rather than 'loop-thru-particles') requires
each particle object to remember the sum of all previously-applied forces as we
compute each successive force; particles do this with a 'force accumulator'
member.

'Force accumulator' members are a much better idea than trying to make forces
into an attribute of particles themselves, because our force-making objects may
be very diverse. Given a force-making object, its easy to say which particles it
affects, but given a particle, it's not always easy to find which force-making
objects will apply to it. For many particle systems we would have to *poll* all
the forces to find which ones contribute to a given particle.  That's OK for
simple systems (e.g. bouncing balls, etc.) but for complex ones such as cloth,
fluids and flocking the number of force-making objects is several times larger
than the number of particles: we don't want to SEARCH for each of the 6 springs
or more that connect to each particle!

***PROBLEM 1****
Particles are supposed to be simple, light-weight, and very, very,very numerous;
we shouldn't design classes/data structs that waste space.
We have a problem, and the Pixar tutorial gives only a clunky solution of
'packing' and 'unpacking' a state vector. BLEAUGH!
What's a better way to tie particles together to make a 'state' vector?

--We DO WANT state vectors to really be a vector: it will let us use other
 people's vector tools, including numerical integration software, (MATLAB,
 Runga-Kutta, verlet, iterative implicit, etc.)

--We DO WANT to avoid slow and tedious copying from particles to state vectors,
from one state vector to another (during integration), from 'next' state vectors
to 'current' state vectors, and copying new state vector results into current
particles.
--We need to update state vectors in a way that resembles double-buffering:
    we use s0 and s0dot to compute s1, then s1 becomes our 'new' s0; repeat.

***SUGGESTION****
A 1-D dyn. allocated array of particle objects *IS ALREADY* a state vector!
    TWO advantages:
    1-as long as the 'particle' objects are all the same size, then an array of
         them in memory forms a nice neat vector, with uniform stride.
    2-you access them only by pointers, so it is easy to 'swap' the contents
    of two particle lists, and thus swap the 'state'; we can 'swap' s0 and s1
    by just changing 2 pointer values, rather than copying all the particle data.

***SUGGESTION***
A CPart class describes the state for one particle.
A dynamically allocated set of CPart objects then forms a state vector.
Note that we can form a DERIVATIVE of a state vector in the same way
    (another dyn. alloc'd array of CPart objects, and define how to add
     and multiply with them).
THUS use data formatted in the same way (a dyn. alloc'd array of CPart objects)
for both the current state vector s0 and its time-derivative s0dot.  This lets
you easily compute Euler integration results:   s1 = s0 + timestep * s0dot.
Later, it lets you compute other solvers easily too
(e.g. midpoint, implicit, RK solvers).

Then, we want to do these things:
        -initialize s0 (the particle system's initial conditions
        -initialize f0 (the list of force-making objects)
        -initialize w0 (the list of constraint-making objects)
        -Make a SINGLE FUNCTION that accepts s0,f0 and w0 to make s0dot.
         it does so analytically: that function describes ALL WE KNOW FOR SURE
         about the trajectory of the particles, or more abstractly, about the
         trajectory of the state s through the (much higher dimensional)
         state space.
        -Make a generic 'solver' function;
            -When 'forward' (Eulerian, explicit, etc.), it computes new state s1
                from timestep, current state s0, and its derivative sdot:
                    s1 = s0 + s0dot*timestep
        -Be ready to make other 'solver' objects(e.g. midpoint,implicit,Verlet).

***IMPORTANT***
Note that DOTFINDER holds the ONLY accurate knowledge we have about the behavior
of our system, of the ideal trajectory of how our state variable will change
over time, and it finds only derivative(s).
Our state changes are all approximations, attempts to find new values that are
on the same trajectory through state space, but advanced by a large timestep.
*/

#endif //end of #ifndef  __CPARTSYS_H__ include guard
